library(readr)
library(tidyverse)
library(ggplot2)
library(Rmisc)
library(DescTools)
library(cowplot)

coda<- read_csv("study4.csv")





coda<-coda%>%
  mutate(age = Q1)%>%
  mutate(male = Q2...2)%>%
  mutate(white = Q3_1)%>%
  mutate(latino = 2-Q4...11)%>%
  mutate(education = Q9)%>%
  mutate(income = Q11)%>%
  mutate(kids=2-kids)%>%
  mutate(bornagain = 2-`born again`)


table(coda$white)

coda$income[coda$income==7]<-NA
coda$education[coda$education==7]<-NA
#coda$white[is.na(coda$white)==TRUE]<-0
coda$male[coda$male==3]<-0
coda$male[coda$male==2]<-0
coda$male[coda$male==4]<-0
coda$pid<-coda$Q2...14
coda$pid[coda$Q2...14 ==1]<-0
coda$pid[coda$Q2...14==2]<-2
coda$pid[coda$Q2...14==3]<-1
coda$pid[coda$Q2...14==4]<-1
coda$pid[coda$Q2...14==5]<-1
coda$daughter<-coda$daughters
coda$daughter[is.na(coda$daughters)==TRUE]<-0
coda$daughter[coda$daughters==1]<-0
coda$daughter[coda$daughters==2]<-1
coda$daughter[coda$daughters==3]<-1
coda$AS1_1[coda$AS1_1==5]<-2
coda$AS1_1[coda$AS1_1==6]<-3
coda$AS1_1[coda$AS1_1==8]<-4
coda$AS1_1[coda$AS1_1==9]<-5
coda$AS1_1[coda$AS1_1==10]<-6
coda$AS1_2[coda$AS1_1==5]<-2
coda$AS1_2[coda$AS1_1==6]<-3
coda$AS1_2[coda$AS1_1==8]<-4
coda$AS1_2[coda$AS1_1==9]<-5
coda$AS1_2[coda$AS1_1==10]<-6
coda$AS1_3[coda$AS1_1==5]<-2
coda$AS1_3[coda$AS1_1==6]<-3
coda$AS1_3[coda$AS1_1==8]<-4
coda$AS1_3[coda$AS1_1==9]<-5
coda$AS1_3[coda$AS1_1==10]<-6
coda$AS1_4[coda$AS1_1==5]<-2
coda$AS1_4[coda$AS1_1==6]<-3
coda$AS1_4[coda$AS1_1==8]<-4
coda$AS1_4[coda$AS1_1==9]<-5
coda$AS1_4[coda$AS1_1==10]<-6
coda$AS1_5[coda$AS1_1==5]<-2
coda$AS1_5[coda$AS1_1==6]<-3
coda$AS1_5[coda$AS1_1==8]<-4
coda$AS1_5[coda$AS1_1==9]<-5
coda$AS1_5[coda$AS1_1==10]<-6
coda$AS1_6[coda$AS1_1==5]<-2
coda$AS1_6[coda$AS1_1==6]<-3
coda$AS1_6[coda$AS1_1==8]<-4
coda$AS1_6[coda$AS1_1==9]<-5
coda$AS1_6[coda$AS1_1==10]<-6


coda$AS2_1[coda$AS2_1==5]<-2
coda$AS2_1[coda$AS2_1==6]<-3
coda$AS2_1[coda$AS2_1==8]<-4
coda$AS2_1[coda$AS2_1==9]<-5
coda$AS2_1[coda$AS2_1==10]<-6

coda$AS2_2[coda$AS2_2==5]<-2
coda$AS2_2[coda$AS2_2==6]<-3
coda$AS2_2[coda$AS2_2==8]<-4
coda$AS2_2[coda$AS2_2==9]<-5
coda$AS2_2[coda$AS2_2==10]<-6

coda$AS2_3[coda$AS2_3==5]<-2
coda$AS2_3[coda$AS2_3==6]<-3
coda$AS2_3[coda$AS2_3==8]<-4
coda$AS2_3[coda$AS2_3==9]<-5
coda$AS2_3[coda$AS2_3==10]<-6

coda$AS2_4[coda$AS2_4==5]<-2
coda$AS2_4[coda$AS2_4==6]<-3
coda$AS2_4[coda$AS2_4==8]<-4
coda$AS2_4[coda$AS2_4==9]<-5
coda$AS2_4[coda$AS2_4==10]<-6

coda$AS2_5[coda$AS2_5==5]<-2
coda$AS2_5[coda$AS2_5==6]<-3
coda$AS2_5[coda$AS2_5==8]<-4
coda$AS2_5[coda$AS2_5==9]<-5
coda$AS2_5[coda$AS2_5==10]<-6

coda$AS3_1[coda$AS3_1==5]<-2
coda$AS3_1[coda$AS3_1==6]<-3
coda$AS3_1[coda$AS3_1==8]<-4
coda$AS3_1[coda$AS3_1==9]<-5
coda$AS3_1[coda$AS3_1==10]<-6

coda$AS3_2[coda$AS3_2==5]<-2
coda$AS3_2[coda$AS3_2==6]<-3
coda$AS3_2[coda$AS3_2==8]<-4
coda$AS3_2[coda$AS3_2==9]<-5
coda$AS3_2[coda$AS3_2==10]<-6

coda$AS3_3[coda$AS3_3==5]<-2
coda$AS3_3[coda$AS3_3==6]<-3
coda$AS3_3[coda$AS3_3==8]<-4
coda$AS3_3[coda$AS3_3==9]<-5
coda$AS3_3[coda$AS3_3==10]<-6

coda$AS3_4[coda$AS3_4==5]<-2
coda$AS3_4[coda$AS3_4==6]<-3
coda$AS3_4[coda$AS3_4==8]<-4
coda$AS3_4[coda$AS3_4==9]<-5
coda$AS3_4[coda$AS3_4==10]<-6

coda$AS3_5[coda$AS3_5==5]<-2
coda$AS3_5[coda$AS3_5==6]<-3
coda$AS3_5[coda$AS3_5==8]<-4
coda$AS3_5[coda$AS3_5==9]<-5
coda$AS3_5[coda$AS3_5==10]<-6

coda$AS3_6[coda$AS3_6==5]<-2
coda$AS3_6[coda$AS3_6==6]<-3
coda$AS3_6[coda$AS3_6==8]<-4
coda$AS3_6[coda$AS3_6==9]<-5
coda$AS3_6[coda$AS3_6==10]<-6

coda$AS4_1[coda$AS4_1==5]<-2
coda$AS4_1[coda$AS4_1==6]<-3
coda$AS4_1[coda$AS4_1==8]<-4
coda$AS4_1[coda$AS4_1==9]<-5
coda$AS4_1[coda$AS4_1==10]<-6

coda$AS4_2[coda$AS4_2==5]<-2
coda$AS4_2[coda$AS4_2==6]<-3
coda$AS4_2[coda$AS4_2==8]<-4
coda$AS4_2[coda$AS4_2==9]<-5
coda$AS4_2[coda$AS4_2==10]<-6

coda$AS4_3[coda$AS4_3==5]<-2
coda$AS4_3[coda$AS4_3==6]<-3
coda$AS4_3[coda$AS4_3==8]<-4
coda$AS4_3[coda$AS4_3==9]<-5
coda$AS4_3[coda$AS4_3==10]<-6

coda$AS4_4[coda$AS4_4==5]<-2
coda$AS4_4[coda$AS4_4==6]<-3
coda$AS4_4[coda$AS4_4==8]<-4
coda$AS4_4[coda$AS4_4==9]<-5
coda$AS4_4[coda$AS4_4==10]<-6

coda$AS4_5[coda$AS4_5==5]<-2
coda$AS4_5[coda$AS4_5==6]<-3
coda$AS4_5[coda$AS4_5==8]<-4
coda$AS4_5[coda$AS4_5==9]<-5
coda$AS4_5[coda$AS4_5==10]<-6


## build sexism measures
coda$hostile_sexism<-coda$AS1_3+coda$AS1_6+coda$AS2_1+
  coda$AS2_2+coda$AS2_3

coda$ben_sexism<-coda$AS1_1+coda$AS1_2+coda$AS1_4+
  coda$AS1_5+coda$AS2_4+coda$AS2_5

coda$ben_men<-coda$AS3_1+coda$AS3_3+coda$AS3_4+
  coda$AS4_1+coda$AS4_3

coda$hostil_men<-coda$AS3_2+coda$AS3_5+coda$AS3_6+
  coda$AS4_2+coda$AS4_4+coda$AS4_5

coda$auth1<-coda$After1-1
coda$auth2<-2-coda$After2
coda$auth3<-coda$After3-1
coda$auth4 <- coda$After4-1
coda$auth5 <-2- coda$After5
coda$auth6<- 2 - coda$After6
coda$auth7<- 2- coda$After7
coda$auth8 <- 2-coda$After8


coda<-coda%>%
  mutate(author = auth1+auth2+auth3+auth4+auth5+auth6+auth7+auth8)%>%
  drop_na(author)


##anova results

RES_mean <- summarySE(coda, measurevar="author", groupvars="AfterBG",
                      na.rm=TRUE, conf.interval=.95)
RES_mean$test<-c("Boy", "Girl", "Neither")
RES_mean


coda.aov<-aov(author~as.factor(AfterBG), data=coda)
summary(coda.aov)
TukeyHSD(coda.aov)


Author_mean_coda <- summarySE(coda, measurevar="author", groupvars="AfterBG",
                              na.rm=TRUE, conf.interval=.95)

Author_mean_coda
Author_mean_coda$test<-c("Boy", "Girl", "Neither")


plot_Author_coda<-ggplot(Author_mean_coda, aes(x=factor(test, level =c('Neither','Boy', 'Girl')), 
                                               y=author, group=1)) +
  geom_line() +
  geom_errorbar(width=.1, aes(ymin=author-ci, 
                              ymax=author+ci), colour="red", 
                data=Author_mean_coda) +
  geom_point(shape=21, size=1, fill="white") +
  xlab("Treatment")+
  ylab("Authoritarianism measure")+ylim(0,5)+ theme_bw()
plot_Author_coda

#measurement results
library(lavaan)
#install.packages("semTools")
library(semTools)
library(irtoys)
library(ltm)
library(mirt)
CFAdata_4<-subset(coda, 
                      auth1 >-1 & auth2 >-1 & auth3 >-1 & auth4+auth5>-1 &
                        auth6>-1 &auth7>-1 &auth8>-1)

vars_4 <- c("auth1", "auth2", "auth3", "auth4",
            "auth5", "auth6", "auth7", "auth8")
irt_data_4 <- subset(CFAdata_4, select = vars_4)

group_4<-as.character(coda$AfterBG )

table(group_4)

group_4 <- dplyr::recode(group_4, "1"="Boy", "2"="Girl", "3"="Neither")

mod_configural_4<-multipleGroup( data=irt_data_4, 
                                 model=1, 
                                 group= group_4 ) 

mod_configural_4
summary(mod_configural_4)
M2(mod_configural_4)

mod_metric_4<-multipleGroup( data=irt_data_4, 
                             model=1, 
                             group= group_4, 
                             invariance=c('slopes'))

mod_metric_4
summary(mod_metric_4)
M2(mod_metric_4)


anova( mod_metric_4, mod_configural_4)




mod_scalar2_4<-multipleGroup( data=irt_data_4, 
                              model=1, 
                              group= group_4, 
                              invariance=c('slopes', 'intercepts', 'free_var','free_means')
)



mod_scalar2_4
summary(mod_scalar2_4)
M2(mod_scalar2_4)

anova( mod_scalar2_4, mod_configural_4)

mod_scalar1_4<-multipleGroup( data=irt_data_4, 
                              model=1, 
                              group= group_4, 
                              invariance=c('slopes', 'intercepts', 'free_var')
)



mod_scalar1_4
summary(mod_scalar1_4)
M2(mod_scalar1_4)

anova( mod_scalar1_4, mod_configural_4)

mod_fullconstrain_4<-multipleGroup( data=irt_data_4, 
                                    model=1, 
                                    group= group_4, 
                                    invariance=c('slopes', 'intercepts')
)



mod_fullconstrain_4
summary(mod_fullconstrain_4)
M2(mod_fullconstrain_4)

anova( mod_fullconstrain_4, mod_configural_4)


anova(mod_metric_4, mod_configural_4) #equal slopes only
anova(mod_scalar2_4, mod_metric_4) #equal intercepts, free variance and mean
anova(mod_scalar1_4, mod_scalar2_4) #fix mean
anova(mod_fullconstrain_4, mod_scalar1_4) #fix variance


anova(mod_fullconstrain_4, mod_scalar1_4, mod_scalar2_4, mod_metric_4, mod_configural_4)


summary(mod_scalar2_4)
coef(mod_scalar2_4, simplify=TRUE)
residuals(mod_scalar2_4)
plot(mod_configural_4)
plot(mod_configural_4, type = 'info')
plot(mod_configural_4, type = 'trace')
plot(mod_configural_4, type = 'trace', which.items = 1:4)
plot(mod_configural_4, type = 'trace', which.items = 5:8)
itemplot(mod_configural_4, 2)
itemplot(mod_configural_4, 2, type = 'RE')




##multinomial logit of what thinking of
##define the baseline
require(foreign)
require(nnet)
require(ggplot2)
require(reshape2)
require(stargazer)
coda$AfterBG<-as.factor(coda$AfterBG)
coda<-coda%>%
  mutate(AfterBG = recode_factor(.x=AfterBG, 
                                 `1`="Boy",
                                 `2`="Girl",
                                 `3`="Neither"))
coda$AfterBG <-relevel(coda$AfterBG, ref = "Neither")


model<- multinom(AfterBG ~ daughter+age+male+white+latino+education+income+
                   bornagain+pid+hostile_sexism+hostil_men+ben_men+ben_sexism, data = coda)
stargazer(model, type="text")

##appendix
table(coda$gender)


sink()